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We review recent work on superlattices in monolayer and bilayer graphene. We high- 
light the role of the quasiparticle chirality in generating new Dirac fermion modes with 
tunable anisotropic velocities in one dimensional (ID) superlattices in both monolayer 
and bilayer graphene. We discuss the structure of the Landau levels and magnetotrans- 
port in such superlattices over a wide range of perpendicular (orbital) magnetic fields. In 
monolayer graphene. we show that an orbital magnetic field can reverse the anisotropy of 
the transport imposed by the superlattice potential, suggesting possible switching-type 
device applications. We also consider topological modes localized at a kink in an elec- 
( mm ) trie field applied perpendicular to bilayer graphene, and show how interactions convert 

these modes into a two-band Luttinger liquid with tunable Luttinger parameters. The 
band structures of electric field superlattices in bilayer graphene (with or without a mag- 
netic field) are shown to arise naturally from a coupled array of such topological modes. 
We briefly review some bandstructure results for 2D superlattices. We conclude with a 
discussion of recent tunneling and transport experiments and point out open issues. 



Keywords: Graphene, Bilayer graphene, Superlattice, Band structure, Transport, Landau 
level, Quantum Hall effect, Luttinger liquid 



1. Introduction 

Graphene is a two-dimensional carbon crystal that exhibits novel physics and trans- 
port properties due to its excitations resembling chiral relativistic massless Dirac 
fermions at low energy.EH^ Its bilayer cousin, Bernal-stacked bilayer graphene (BLG) 
has also garnered much interest due to the possible novel broken symmetry states 
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it could potentially exhibit in the presence of interactions that destabilize the 
quadratic band touching point present in a minimal tight-binding modelJ^Hl Both 
graphene and bilayer graphene are also widely regarded as viable materials for de- 
veloping new types of device applications due to the chiral nature of their low energy 
excitations. This is largely due to the possibility of opening tunable band gaps - by 
engineering a relative potential difference on each sublatticeP^HIH or by strain en- 
gineering) 16 ! 17 ! in monolayer graphene or by applying an electric field perpendicular 
to the layers in bilayer graphene.'^^ni 

In this review, we focus on new physics that develops in the presence of slow 
spatial potential variations in both monolayer and bilayer graphene. While a spa- 
tially varying chemical potential could result in p-n junctions in either case,^2H 
a perpendicular electric field that changes direction as a function of position, spe- 
cific to bilayer graphene, leads to midgap domain wall states (or 'kink' states) that 
have topological character J22HHI For an isolated "nanowire" , electron interactions 
drive this ID system into a tunable two-band Luttinger liquid.^ These tunable 
nanowires have also been shown to act as controllable "electronic highways" In 
both monolayer and bilayer graphene, profiles with periodic potential variations 
form superlattices that can lead to dramatic modifications in the bulk band struc- 
ture. Such superlattice potentials have been shown to induce new Dirac modes at 
zero and finite energy with tunable velocities and transport propertiesJ^SHiSli!] 

We review the physics of such superlattices as well as the effect of a perpendic- 
ular orbital magnetic field for three experimentally accessible field regimes: weak, 
moderate, and strong fields.^ A weak orbital magnetic field essentially acts as a 
'probe' of the Dirac modes, while a strong magnetic field overwhelms the effect of the 
superlattice potential leading to quantum Hall physics indistinguishable from that 
of pure graphene. At moderate magnetic fields, however, we find dispersing Landau 
levels and interesting field tuning of transport properties. In addition, we discuss the 
effect of a magnetic field on the so-called topological 'kink' states that form at the 
interface that separates two region with opposite interlayer bias.^JEZl Throughout 
this article we discuss recent developments, ongoing experimental efforts, and open 
issues. 

2. Superlattices (SLs) in monolayer graphene (MLG) 
2.1. Bandstructure of ID superlattices 

For pristine MLG, ignoring spin, the low energy Hamiltonian is given by a 2 x 2 ma- 
trix at each valley, Ho = Vf(sp x o~ x — p y (7y), where pseudospin a z — ±1 labels the two 
trigonal sublattices, while the two (decoupled) valleys at ±K = ±4irx/\/3a are la- 
belled by s = ±l. Here, Vf — 3ta/2 is the isotropic Fermi velocity, with a — 1.42 Aand 
t = 3eV being the nearest neighbor carbon-carbon distance and transfer integral re- 
spectively, p is the momentum measured from K. (We set h = l for convenience.) 
The quasiparticles in the vicinity of each valley then behave as massless linearly 
dispersing Dirac fermions, with an energy dispersion U/|p|. 
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We focus here on the effect of a smooth SL potential, for which the period of 
the SL is significantly larger than the interatomic distance. This means we can 
safely ignore the intervalley scattering of electrons which involves large momentum 
transfer. We therefore use the above low energy Hamiltonian and focus on the 
electronic properties of SLs near a single valley. To this end, a ID SL potential 
can be modelled as Hsl = U(y)I, where U(y) = U(y + A) with A being the SL 
period, and / is the identity matrix in the pseudospin (sublattice) space. To gain a 
qualitative understanding of the nontrivial phenomena arising from such SLs, such 
as the anisotropic Fermi velocity renormalization or SL induced band gaps, it is not 
necessary to assume any specific form for U(y). 

The problem of finding the energy spectrum of H = Hq + Hsl has been ex- 
tensively studied. It was noticed by Park et. alJ^ that the chiral nature of Dirac 
fermions in graphene leads to the anisotropic renormalization of Fermi velocity near 
Dirac point. Surprisingly, the Fermi velocity is not renormalized in the SL direc- 
tion, but is suppressed perpendicular to the modulation direction, a counterintuitive 
effect that is deeply rooted in the chiral nature of the Dirac fermions in MLG. 



2.1.1. Weak SL Potential 



For a weak SL potential, the energy spectrum near Dirac point and Fermi velocity 
renormalization can be well understood from perturbation theory. By expanding the 
Hamiltonian H in the chiral basis, |ps) = ^(1, se~ l6k ) T , where cos 9 P — p x /\p\ and 
s = ± denotes electron and hole states, the kinetic energy part Ho can be brought 

with 

), 



ss'e PiP+ 



into diagonal form, while the matrix elements between |ps) and |p + nG, 
G = (0, 2tt/\) as the reciprocal lattice vector, is given by C/ (" G ) (\ 
where pp+n G = 9 p — p + n G- Therefore, for states with momenta parallel to the SL 
direction, p = (0,p v ), we can show that the full Hamiltonian matrix will consists of 
two decoupled blocks, 

/-• \ 

(p-G) U(G) U(2G) 
U*(G) e (p) U(G) 
U*(2G) U*(G)e h (p + G) 

CD 

e e (p + G) [/(G) U(2G) 
U*(G) e h (p) U(G) 
U*(2G) U*(G)e h (p-G) 

\ ••/ 

where e e / l (k) = ±Up|k| are the electron (hole) energies. Applying the second order 
perturbation theory, the energy correction at momentum k is given by 



E 



\U{nG)\' 



(2) 
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where the sum is carried out in the same block and e(p) is understood as the 
corresponding electron or hole energies. This term is zero because of the linearity 
of the spectrum, i.e., e(p) — e(p — nG) = e(p + nG) — e(p). This means, along the 
direction of the superlattice, the Fermi velocity is not renormalized. Therefore, the 
absence of Fermi velocity renormalization in the SL direction is a consequence of 
the chiral nature of Dirac electrons and the linearity of the spectrum. 

At the MBZ boundary, p = (0, ±7r/A), the two blocks become exactly identical, 
which means the energy will be doubly degenerate and energy spectrum is gapless 
at this point. This result is exact and independent of perturbation theory, and this 
band touching point will always be present. 

Once the momentum p is no longer parallel to the reciprocal lattice vector 
G, this nice decoupling will break down and Fermi velocity in the corresponding 
direction, vp = v(p) • p, will inevitably become renormalized. Within second order 
perturbation approximation, the renormalization with respect to pristine MLG is^ 

v v~ v f _ 2\U(nG)\ 2 2 



Vf n 2 v 2 f \G\ 2 

1 n^O /' 



sin^ p . G , (3) 

where # Pi g is the angle between p and G. Since the right hand side of Eq. ^ 
is always negative, except for p || G, the Fermi velocity is decreased from pristine 
MLG value. In contrast, for artificial electrons with linear dispersion but no chirality 
in a weak ID SL, the second order perturbation result for the Fermi velocity is given 

by 

vp-vf = _ V 2\U(nG)\ 2 

which is isotropically decreased and is independent of the direction of p. Therefore, 
the anisotropic renormalization of the Fermi velocity near Dirac cone is truly a 
signature of chiral low energy excitations in MLG. The left Figure [l] shows the 
dispersion for a weak ID SL potential, and corroborates the above results. 

For a square barrier SL (Figj2|, the energy spectrum can be found exactly.^^l 
By making use of Bloch theorem and matching boundary condition, it can be shown 
that the energy spectrum can be obtained from the following transcendental equa- 
tion, 

cosp x = cos(X w l w ) cos(A fc / h ) - Qsm(X w l w )sin(\ b l b ). (5) 
Here, we have used the following notation: 

Syj — S ~\~ 111 fa j Si) — - S lllnjy 



UoX _ W b ,, 

v f A 



\ w = (el-pl) 1/2 , X b =(st-pl) 1/2 and Q = E <° e » & . ( 6 ) 

A w A b 

For a symmetric SL, W b = W w , or equivalently lb = l w = |, Eq. (I5| reduces to 

Am A^ _ . X w , X b ,_v 

cos p r = cos — cos Q sm — sm — , 7 

yx 2 2 2 2 
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Fig. 1. Left: Energy spectrum for a single anisotropic Dirac cone. [Reprinted by permission from 
Macmillan Publishers Ltd: Nature Physics 4, 213, (2008).] Notice the robust band crossing and 
nontrivial minigap opening at the MBZ boundary. Right: Energy spectrum showing five anisotropic 
Dirac cones. [Reprinted with permission from RefGS] Copyright (2010) American Physical Society], 



2U 




Fig. 2. Step SL potential, where Vfo is the SL potential strength, A is the SL period, and Ww(Wb) 
is the width of potential well (barrier). 



where e w = e+u/2 and £{, = e—u/2. For this symmetric case, the energy spectrum is 
particle-hole symmetric, which means Eq. ^ is invariant under the transformation 
£ — » — E. 

To obtain the behavior near K point, we can expand Eq. Q in small e and p x . 
The result is 

From this, we can see that the low energy spectrum is indeed described by an 
anisotropic Dirac cone, with v y = v/ and 



(9) 
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2.1.2. Strong SL potential 

Since the Fermi velocity perpendicular to the SL direction can be significantly renor- 
malized and even brought to zero for a broad region in momentum space, the energy 
spectrum becomes dispersionless in this direction and the electrons can be colli- 
mated in the SL direction.^ Moreover, extra Dirac points can be generated in the 
energy spectrum for an even stronger ID SL, which are shown in the right panel of 

Fig |^ I30I33I36I38| 

To determine the condition for the emergence of extra Dirac points and also 
their locations, we consider a symmetric SL and assume p y = and e = in Eq. 
0. Then Eq. Q reduces to 

2 /A t0 \ u 2 /4 + p 2 2 (X W \ 
1 = cos' — + -^j- A 1 sin' — , 10 

which can be solved by either w 2 /4 + p 2 = w 2 /4 — p x or sin 2 (A u) /2) = 0. The former 
condition gives p x = 0, which is just the original Dirac point. The latter condition 
leads to A^/2 = jit with j being a nonzero integer. Then, the position of new Dirac 
points are subsequently found at 



p x = ± ] j^-4jW. (11) 

For asymmetric SL, the energy spectrum is no longer particle hole symmetric and 
extra Dirac points will appear with nonzero energies. 

At the induced Dirac points, the Fermi velocities behave differently from the 
original Dirac point.^ To see this, we can expand Eq. Q in e up to second order 
and obtain 

1/2 



e± =± 



4|a 2 | 2 [fc 2 sin 2 (a/2) + a 2 sin 2 (fc y /2)] 



(12) 



k^asina + a 2 w 4 /16 - 2k 2 x u 2 sin 2 (a/2). 

with a = [u 2 /4 — fc 2 ] 1 / 2 . Then, at the jth extra Dirac point, the Fermi velocities 
along x and y directions are given by 

u 2 /4-4j 2 ^ 2 16^ 2 j 2 cos(fc a /2) 
v x = v f , v y = v f . (13) 

In contrast, at the original Dirac point, v x — Avf sin(u/4)/u and v y = Vt. 
2.2. Zero field transport 

By assuming a constant relaxation time at the Fermi energy r(Ep) — Tp 1 the dc 
conductivity for MLG SL can be calculated by^ 

a u (E F ) = ^ ~ (14) 

n.k 

where v n i — (nk\vi\nk) is the average velocity in the i-th direction for n-th energy 
band, / n k = l/{exp[f3(E n \ l — Ep)] + 1} is the Fcrmi-Dirac distribution with f3 = 
l/k B T. 
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Fig. 3. (a), (b) a yy and a xx as a function of Fermi energy for an SL with u = 6ir, for Z;, = 0.5 
(red curve) and ij, = 0.4 (dashed blue curve) respectively. (c),(d) a yy and (T xx as a function of 
Fermi energy with Zj, = 0.5 for different SL potential strength. [Figures reprinted with permission 
from Ref.l2£J. Copyright (2010) American Physical Society.] 



The results for the various conductivities as a function of Fermi energy are 
shown in Fig. [3j Fig. [3j[a) and[3^b) show a yy and cr xx , respectively, for an SL with 
u — 6ir and /3 = Vf/ksTX = 20. Here, red and dashed blue curve correspond to 
symmetric (lb = 0.5) and asymmetric (lb = 0.4) SL, respectively. From Fig. [3]ja), 
notice that a yy is oscillating when the Fermi energy is below the barrier but increases 
on the average almost linearly when the Fermi energy is above the barrier. On the 
other hand, a xx always increases on the average with the Fermi energy. Both a yy 
and a xx show oscillating behavior and are symmetric (asymmetric) for symmetric 
(asymmetric) SLs. For <J yy , however, there is a dip at the crossing energies of those 
mini bands. From Eq. |5]), it can be shown that crossing energies occur at e = nn 
at p x = 0. 

Fig.J^c) and|3^d) show a yy and a xx respectively for a symmetric SL with differ- 
ent SL potential strength. Notice that, at low energies, a xx is smaller than its value 
in the absence of an SL. Upon increasing SL potential strength, a xx also increases 
due to the appearance of extra Dirac points. Also, we can notice that a xx is always 
smaller than a yy at low energies, since v y > v x near the Dirac point. 

Now let us consider conductivities for symmetric SLs at zero temperature and 
charge neutrality (i.e.,T — and Ep = 0) with only one Dirac point in the spectrum. 
As the SL potential is not very strong, the low energy Hamiltonian can be written 
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as H = v x k x a x + v y k v a y , with anisotropic Fermi velocities. It can be shown that 
the conductivity along and perpendicular to the SL direction is given by 

ayy{ E F = 0) = -a = <?o ^ , a xx (E F = 0) = -a = ^o^^, 

(15) 

where <jq is the universal conductivity of an isotropic Dirac cone. The value of <jq 
depends on the order of different limits taken, such as vanishing temperature and 
Fermi energy. However, the form of the result for an anisotropic Dirac cone does 
not depends how <tq is calculated. 

When there are extra Dirac points in the spectrum, by assuming their indepen- 



dence and using Eq. (|15|), the conductivities now are 

|sin(u/4)| , , («/4) 2 -iV 



a vy (E F = 0)=a a \ — + 2^ , 



o xx {Ef = 0) = a I " /4 + 2 + 7 /A l n 22 }' ( 16 ) 
' |sm(M/4)| ^—^ [u/iy — 2 ^ 



where j max = Int[u/47r] counts pairs of extra Dirac points. From Eq. ( 16 ) we can see, 
every time a new pair of Dirac points is generated in the spectrum with u — 4mr and 
n an integer, the conductivity parallel to the SL will shows a dip, while conductivity 
in the perpendicular direction will diverge. 

These results can be confirmed by numerical calculation, using Landauer- 
Biittiker formalism.^1 Fig. |4] shows the corresponding numerical results. In the left, 



the conductivity parallel to the SL direction (top panel) agrees with Eq. ( 15 ) when 
U$ is smaller than the critical value where extra Dirac points are generated. The 
corresponding Fano factor is 1/3 (bottom panel), which agrees with pseuodiffusivc 
character of transport. Once Uq exceeds the critical value, new Dirac points will 
emerge and provide new transmission channels in the SL direction. We can see that 



the conductivity qualitatively agrees with Eq. (16). Since Eq. (16) is based on the 
assumption that all the Dirac points are independent and have linear dispersion, 
which is valid in a very small energy region, it is not surprising to see the numerical 
results depends on both SL potential strength and SL period. Also, the Fano factor 
is larger than 1/3, indicating that the transport is no longer pseudodiffusive. 

In contrast, conductivity perpendicular to the SL direction is shown in the right 
figure of Fig. [4] Again, when SL potential strength is smaller than the critical 



value, the conductivity is well described by the simple picture of Eq. ( 15 ). When U 



approaches the critical value, the conductivity shows a peak. For even stronger SL 



potential, the numerical result agrees with Eq. ( 16 1 quite well, which may suggest 



that the approximations adopted for Eq. ( 16 1 are appropriate in the perpendicular 



direction. Also, the Fano factor is 1/3 for almost all SL potential strengths, except 
for those critical values where new Dirac points are generated. 



March 6, 2013 11:34 WSPC/INSTRUCTION FILE GrapheneUMPB 



Graphene: Kinks, Superlattices, Landau Levels, and Magnetotransport 9 



13 4 



2 - 



0.4 - 



0.3 



1 1 1 1 

(7|| 


i i 


- A =24 a 




~~ A =36 a 

1 | 1 -f 


/ # — 

- — 1 










2;r 



Uq a/ hv i 



3ji 



40 

-~ 30 
^ 20 

H 

10 



0.3 



0.1 




UoX/hvi 



Fig. 4. Conductivity parallel (left) and perpendicular (right) to the SL direction and the corre- 
sponding Fano factors. Here, A is the SL period and a = 1.42 Ais the lattice constant of graphene. 
The black dotted line corresponds to conductivity calculated from Eq. JlM (top panel) and Fano 
factor F = 1/3 (bottom panel). [Reprinted with permission from Ref.QSl. Copyright (2011) Amer- 
ican Physical Society.] 



2.3. Landau levels 



In a uniform perpendicular magnetic field, the eigenenergy for pristine MLG in the 
absence of SL is e n = sgn(n)yjnju; c , where ui c = ^/2vp / £b, with £b = 1/yeB. For 
s = +1 (i.e., at valley K), the eigenfunctions are given by 

Akx 



<Pn,k,+ (x,V) 



1p\n\,k(y) 



V2L \sgn(n)i)\ nl _ ltk (y) 



(17) 



where L and k are the system length and electron momentum deviation from K, 
both along the x-direction, while for n = 0, 

Akx 



o 



Here, ip n ,k(y) is the n-th eigenstate of a (shifted) ID harmonic oscillator, 

2" 



1 



^AV) = 7r , ex P 



i / y - 2/0 

2 



y-yp 
£ B 



(18) 



(19) 



centered at yo — k£ B , and H n are Hermite polynomials. For s — —1 (i.e., at — K), 
the eigenfunctions are given by 4> n .k-{ x ,y) = ^i (J y4>n,k,+ {x,y). The full low energy 
LLs of MLG are thus 4> n ,k,±{ x -> y)e ±iK: ° x . 

We now turn to the effect of a periodic ID chemical potential modulation V(y), 
with period A>a, on these Landau levels at low energy. Recent work has shown 
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these results to be generally consistent with solving the Harper equation using 
the full tight-binding model.^ The set of eigenfunctions <j) n ,k,s{x, y)e lsKxX , with 
s = ±1, form a convenient basis to study the SL Hamiltonian in a magnetic field. 
(This basis choice is different from the one used by Park, et alj^ and allows us 
to numerically access a wide range of magnetic fields.^ In the weak field regime, 
our results are consistent with Ref.^H) Due to momentum conservation along the 
x-direction, the SL Hamiltonian is diagonal in k. Further, for ) » a, intervalley 
scattering is strongly suppressed. We will therefore assume that the two valleys 
stay completely decoupled. (We focus below on valley K with s = +1; we expect 
identical physics around valley — K.) With this approximation, the only effect of 
the SL potential is, thus, to induce Landau level mixing. 

To proceed, we need to choose a concrete form for the SL potential. For sim- 
plicity, we set V(y) = ^ cos (^p) , although our results can be easily generalized to 
other (e.g., step- like) SL potentials by including multiple Fourier components. We 
can then expand the Hamiltonian in the above basis, retaining up to 3000 Landau 
levels, and diagonalize it to obtain the spectrum of the ID SL in a magnetic field. 

In order to study the effect of the magnetic field on the ID SL in graphene, with 
U = U\/2nvf ~ 0(1), it is useful to consider three regimes for the magnetic field. 

(i) Weak field: This regime corresponds to having fvjj c -C U, where the Landau 
level spacing is much smaller than the SL amplitude, so that 2it£b/X 3> 1. In this 
regime, the magnetic field may be viewed as effectively 'probing' the zero field SL 
excitations. 

(ii) Intermediate field: In this regime, Hu) c ~ U, which means 2tt£b/X ~ L so 
that the SL potential and the magnetic field have to be treated on equal footing. 

(iii) Strong field: Here, huj c U or, equivalently, 2tt£b/^ <C 1. In this regime, 
the SL potential only weakly perturbs the Landau levels of pristine graphene. 

Fig. [5] shows the spectrum of the graphene SL in different field regimes for SL 
strengths U=2irv / / 'A (or U = 1) and 6iTVf/X (or U = 3). This allows us to contrast 
the behaviour of the spectrum of the SL in a magnetic field without or with extra 
Dirac points being present at zero field, and to explore consequences for quantum 
Hall physics and transport. 

2.3.1. Weak field regime 

When the magnetic field is weak, £g = 2A (top panels in Fig. [5J, we find that the 
energy spectrum barely depends on the value of k, or equivalently, j/o- This is due 
to the fact that when magnetic length £b is larger than the SL period A, the matrix 
elements of the Hamiltonian do not depend on the center of the LL wavefunctions, 
which yields flat bands. Equivalently, in this regime, the magnetic field may be 
viewed as effectively 'probing' the structure of the zero field SL dispersion leading 
to Landau levels which depend on the nature of the Dirac spectrum at low energy. 

For U = 1, the low energy spectrum of the SL contains a single anisotropic 
Dirac point at zero energy. For an anisotropic Dirac cone described by an effective 
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Fig. 5. (Color online) Landau levels of monolayer graphene SL for different (dimensionless) SL 
strengths U, and magnetic fields B. The spectrum is shown for weak field (£g = 2A, top panels) 
and intermediate field (Is = 0.2A, bottom panels). Left panels (a,b) correspond to U = 1 which 
supports a single anisotropic zero energy massless Dirac fermion. Right panels (c,d) correspond to 
U = 3 which supports three zero energy massless Dirac fermions with anisotropic velocities - the 
weak field zero energy LL thus has three times as many states for U = 3 as it does for £7 = 1, while 
the n = ±1,±2 levels have degeneracy splitting in weak field due to the Dirac fermions having 
two different mean velocities. For !j «A (not shown), the LLs closely resemble that of pristine 
graphene. See text for a detailed discussion of the Landau level structure. 



Hamiltonian H — v x k x a x +v y k y a y , the LLs are given by e n — sgn[n)yj2\n\v x v y /iB- 
Since the SL renormalizes v x < Vf, but leaves v y — Vf, the Landau levels at weak 
field resemble those of pristine graphene, but with a renormalized lower effective 
velocity ^Jv x v y < Vf. 

For U = 3, the low energy spectrum of the SL contains three anisotropic Dirac 
points at zero energy, so that the zero energy Landau level has three times the 
degeneracy of the case with {7 = 1. Further, the Dirac cone centred at K has a 
slightly different average velocity ^Jv x v y compared with the two cones which are 
symmetrically split off from K along ±x. This degeneracy breaking results in the 
Landau levels at nonzero energy becoming weakly split, as is most clearly seen 
for the first two excited Landau levels (at positive or negative energy, i.e., with 
n = ±1, ±2). We have numerically determined v x and v y for each of the three Dirac 
points and found good agreement between the energy levels obtained on this basis 
of having Dirac fermions with two different average velocities, and that obtained 
directly numerically. 
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At higher energies, E/u c > 2 for U = 1 or E/u> c > 1 for {7 = 3, the spectrum 
begins to deviate from this simple behavior expected for a linear Dirac spectrum. 
This deviation results from curvature in the dispersion, which appears upon going 
beyond the linearized approximation. 

2.3.2. Intermediate field regime 

At intermediate fields, for £b = 0.2A, the spectrum at low energy is most simply 
understood as arising from the SL potential inducing a strong dispersion to the 
Landau levels. In simple terms, if we assume that the state labelled by momentum 
k, or equivalently position yo, have an energy which is modulated by the SL poten- 
tial, we expect a periodic modulation of this energy with period A and amplitude 
proportional to the SL amplitude U. The behaviour of the low energy Landau lev- 
els, n = 0, ±1,±2, as seen from the lower panels in Fig[5j is consistent with this 
scenario, with the modulation following the cos(27ry/A) form of the SL potential 
and the modulation for U — 3 being roughly thrice as strong as the modulation 
for U — 1. We can also see that the low energy Landau levels when U — 3 overlap 
with each other. This will have nontrivial effect on the dc conductivity, as shown in 
the following subsection. For higher energy Landau levels, the energy spectrum still 
has a periodic modulation but no longer retains the simple form of cosine function. 
This is due to the fact that as the energy gets higher, the distribution of Landau 
levels becomes more dense and the energy difference between two adjacent levels 
is now comparable to the matrix element of SL potential. Therefore, a simple first 
order perturbation correction is not enough to account for the dispersion and second 
order perturbation from adjacent levels must be taken in account, which causes the 
Landau level to lose its simple cosine form. 

2.3.3. High field regime 

For very strong magnetic field, the Landau level structure of pristine graphene is 
recovered. Here, only one zero energy level exists at the Dirac point, and other en- 
ergy levels follow the square root relation. This is simply because in such a strong 
magnetic field, the SL is just a perturbation and can only give rise to a small modu- 
lation of the LLs following our argument at intermediate field. From a perturbative 
point of view, the energy corrections up to first order to the LL energies are given 



which gives a sinusoidal dependence on the center position of LL wavefunctions. 
Thus, even in a strong magnetic field, the energy spectrum is not dispersionless but 
has a spatial modulation following the SL; however the ratio of the amplitude of 
this modulation to the Landau level spacing, U /cj c , is extremely small in the high 
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Fig. 6. (Color online) Diagonal dc conductivities of monolayer graphene SL for different (dimen- 
sionless) SL strengths U, and magnetic fields B. The conductivity is shown for weak field (£g = 2A, 
top panels) and intermediate field (£g = 0.2A, bottom panels). Left panels (a,b) correspond to 
U = 1, and right panels (c,d) correspond to U = 3. The conductivities show strong anisotropy 
when magnetic field strength is tuned - for weak field (a,c), a yy is larger than cr xx , which is a 
consequence of the Fermi velocity renormalization in the absence of magnetic field; for moderate 
field (b,d), the anisotropy is reversed, since v x acquires intra-LL contributions, as explained in 
the text. For £g A (not shown), result for pristine graphene is recovered and the transport is 
isotropic in both directions. 



field regime. This dispersion, though small, can give rise to interesting magnetoresis- 
tance oscillation known as Weiss oscillation, on top of the usual Shubnikov-de Hass 
oscillation.^! It was shown that, compared to two-dimensional electron gas with 
parabolic dispersion relation, Weiss oscillation in graphene SL is more pronounced 
and is more robust against temperature damping in small field region. This is a 
consequence of the different Fermi velocities of Dirac and normal electrons at same 
chemical potential.^ 



2.4. Magnetotransport 

Once we have the eigenvalues and eigenf unctions for the superlattice in a perpen- 
dicular magnetic field, both ac and dc conductivities can be calculated directly by 
Kubo formula, 

M = e2 1 f x d y f( E °) - /Cgg) («fchl/^)(/3fcN«fc) , v 

aiAUJ} hjr\e 2 J Va E a — Ep E a -Ep-w-iT ■ { ' 
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(a) (c) 




Fig. 7. (Color online) The dc Hall conductivity of monolayer graphene SL for different (dimen- 
sionless) SL strengths U, and magnetic fields B. The conductivity is shown for weak field (£g = 2A, 
top panels) and intermediate field (£b = 0.2A, bottom panels). Left panels (a,b) correspond to 
U = 1, and right panels (c,d) correspond to U = 3. For weak field (a,c), the Hall conductivity 
shows well-defined plateaus, as a consequence of nearly flat energy bands. For intermediate field 
(b,d), the energy bands become dispersive and the Hall conductivity no longer shows step-like 
structure. However, for weak SL (b), the energy bands are not fully overlapped, Hall conductivity 
still shows small plateaus when chemical potential falls between two bands, and the value of a xy 
changes by one between adjacent steps, as expected from Dirac physics. For £g <C A (not shown), 
result for pristine graphene is recovered and Hall conductivity is constant between adjacent LLs 
and changes by one when chemical potential crosses an LL. 



Here, we have set T = 10~ 3 x 2irvf/X as the Landau level broadening, E a (yo) 
and \ak) are the a-th eigenvalue and the corresponding eigenstate of the system 
which can be expanded in the basis of \nk) 1 where 4> n (k,y) — (y\nk) is the LL 
wavefunctions for pristine graphene. Vi is the velocity operator in i-direction and 
t>i = vpo'i-, where di is the Pauli matrix. Note that (ak\v y \ak) — is always true 
for any state. 

Results for dc diagonal conductivities as function of chemical potential fj, are 
shown in Fig. [6j This can be done by setting the frequency w to zero in Eq. (21 ), 
and only the real part of the conductivity tensor is nonzero. In weak magnetic 
field, the conductivities show strong anisotropy, with a yy larger than a xxi which 
is a consequence of the Fermi velocity renormalization in the absence of magnetic 
field (see Fig. [6] (a) and (c)). Since {ak\v y \ak) = and (ak\v x \ak) ~ because 
of the flat band structure, the major contribution to the diagonal conductivities 
comes from off-diagonal matrix elements, (ak\vi\(3k) with a ^ /3. Numerically, we 
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have observed that matrix elements of v y is always larger than those of v x , which 
gives rise to the anisotropy in the weak field. In intermediate magnetic field, con- 
ductivities still show anisotropy, but with <j xx significantly larger than a vy (see 
Fig. [6] (b) and (d)). This is because v x has acquired diagonal matrix element, 
(ak\v x \ak) = dE a (y Q = kl 2 B )/dk since the energy spectrum is dispersive, 
while v v still lacks this contribution. Notice the positions of the conductivity peaks 
of (jyy exactly correspond to the minimum and maximum of the energy band, where 
the density of states is the largest. For a xx , however, the conductivity is minimum 
at the band edge, since the average of the velocity operator, (v x ), is zero. Therefore, 
the intra-LL contribution to a xx is the smallest at the band edge. For weak SL 
potential, a xx can drop to zero when there is no overlapping LLs, while in a strong 
SL, a xx always show dispersive transport property. In strong magnetic field (not 
shown here), where the Landau levels structure of pristine graphene is recovered, 
the conductivities become isotropic (see, for example, Refill). In this case, the SL 
is merely a perturbation to the magnetic field and thus should have minor effect on 
determining the magnetotransport properties. 

The dc Hall conductivity is shown in Fig. [7j For weak magnetic field (Fig. [7] 
(a) and (c)), the Hall conductivity shows well-defined plateaus, as a consequence of 
nearly fiat energy bands. The values of Hall conductivity around Dirac points are 
±l/2(e 2 /h) in weak SL (U = 1) and ±3/2(e 2 /h) in strong SL (U = 3). This result 
resembles the anomalous half integer quantum Hall effect in pristine graphene and 
the Hall conductivity triples due to the existence of three Dirac points in a strong 
SL. Moving away from the Dirac point, we can observe quantum Hall plateaus with 
higher conductivities, and the value increases by 1 each time the chemical potential 
crosses an LL. For intermediate magnetic field (Fig. [7] (b) and (d)), there is no 
longer well defined plateaus due to the dispersive energy spectrum. However, for 
weak SL, the LLs are not overlapped with each other. If chemical potential falls 
between two LLs, a small plateau can still show up, with the value expected from 
Dirac physics. When magnetic field becomes strong enough as the LL structure for 
pristine graphene is restored, Hall conductivity will show anomalous half integer 
quantum Hall plateaus. 

Fig. [8] shows the ac conductivities of graphene SLs in an intermediate magnetic 
field. For weak and strong magnetic fields, the results resemble those of pristine 
graphene,^ since in both cases the LLs are nearly flat and the real part of the con- 
ductivities show strong peaks when photon frequencies exactly correspond to the 
energy differences between two LLs. In an intermediate magnetic field, the result is 
complicated by the dispersion of LLs. At low frequencies, there can be optical tran- 
sitions in a range of photon energies, and the real part of diagonal conductivities 
is maximum at the band edge where the DOS is also maximum. At high frequen- 
cies, the LLs become less dispersive and peaks will show up. These results can be 
linked with graphene's unusual magneto-optical properties, for example, giant Fara- 
day rotation.! 50 ! 51 ! While the anisotropy in the diagonal conductivities can lead to 
anisotropic rotation angles for incident waves with different polarization plane, this 



March 6, 2013 11:34 WSPC/INSTRUCTION FILE GrapheneUMPB 



16 M. Killi, S. Wu, and A. Paramekanti 




Fig. 8. (Color online) The ac conductivity of monolayer graphene SL for different (dimensionless) 
SL strengths U, in intermediate magnetic field, £g = 0.2A, with fi = 0.2ui c . Left panels (a,b) 
correspond to U = 1, and right panels c,d) correspond to U = 3. 

effect is actually quite small and hard to observe experimentally. 



2.5. Bandstructure of 2D superlattices 

In 2D SLs, the Fermi velocity near the Dirac point is anisotropically renormalized 
along every direction. Due to the chiral nature of low energy excitation, there are 
still energy band crossing at the MBZ boundary. 

Fig. [9] shows a 2D rectangular SL with muffin tin type SL potential with pe- 
riod L x and L y in x and y directions, and the corresponding energy spectrum. In 
contrast to ID SL where Fermi velocity parallel to the SL direction is not affected, 
Fermi velocity in a rectangular SL is renormalized in every direction. This can be 
clearly demonstrated by second order perturbation, assuming a weak SL potential 
strength,^ 

1 



27T 2 [/ 2 V 2 _ 



— Jsm^ k)G , 



(22) 



where U^d is the 2D SL potential strength in a circular region of diameter d, G = 
(2nm/L x , 2Trn/L y ) is the reciprocal lattice vector with m and n integer, and J\{x) 
is the Bessel function. Since G can be along any direction, compared to ID SL 
where G is always along the SL direction, we can see that the Fermi velocity is 
renormalized in every direction. 
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Fig. 9. (Color online) 2D rectangular muffin-tin type SL potential leads to anisotropically renor- 
malized Dirac cone (left) but no minigap at the MBZ boundary (right). [Reprinted by permission 
from Macmillan Publishers Ltd: Nature Physics 4, 213, (2008).] 



The energy spectrum of a 2D rectangular SL also has band crossing points 
in the middle of an MBZ boundary edge, similar to ID SL. In addition to these 
crossing points, at the four corners of the MBZ, energy gap also closes. When 
similar calculation is carried out for an artificial non-chiral electron, these band 
crossing points disappear, which truly suggests that they are the consequence of 
the chirality of low energy excitations in MLG. 

Even though band crossing points appear at the MBZ boundary in both ID and 
2D rectangular SLs, the density of states does not vanish at the crossing energy and 
the newly generated massless Dirac fermions are obscured by other states. However, 
for triangular SLs, there exists an energy window where the only available states 
come from the newly generated massless Dirac fermions.^ Fig. 10 shows a trian- 
gular SL with muffin-tin type SL potential and its corresponding energy spectrum. 
Again, the Fermi velocity is anisotropically renormalized in every direction. The 
gap between the first and second conduction bands vanishes in the middle of the 
MBZ boundary edges, and the density of states also vanishes linearly here. This 
result will have significant impact on the experiments explained later J 3 -^ 37 * 43 * 52 ! 

When graphene is expitaxially grown on a substrate (i.e., SiC, hexagonal 
boron nitride (hBN), transition metal surfaces, etc.), the lattice mismatch between 
graphene and the substrate and also their relative orientation can lead to a 2D SL 
with large period. Therefore, theoretical results can be tested on these structures. 
However, previous results are based on an effective Hamiltonian approach which 
assumes that external potential does not break sublattice symmetry. When such a 
symmetry breaking effect is taken into account, most of the earlier results will be 
modified. For example, a gap should open up at Dirac point and minigap should 
appear at the MBZ where bands are backfolded. Pletikosic et. al.^ have observed 
a minigap in graphene expitaxially grown on Ir(lll) surface, which is due to Moire 
patterned periodic potential. They could not determine whether the Dirac point is 
gapped because graphene on Ir(lll) is slightly p-doped. On the other hand, Rusponi 
et. al.^ showed that, in the presence of sublattice symmetry breaking SL potential, 
the Dirac point remains intact and, remarkably, the Fermi velocities are anisotrop- 
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Fig. 10. (Color online) 2D triangular muffin-tin type SL potential (left) gives rise to finite 
energy massless Dirac fermions (right). [Reprinted with permission from RefM Copyright (2008) 
American Physical Society.] 

ically renormalized and the energy spectrum becomes trigonally warped. This is 
consistent with the theory of Ortix et. al.,^1 where it was demonstrated, incommen- 
surate Moire patterned superstructure preserves the Dirac cone in a renormalized 
form, with threefold global symmetry due to a substrate-induced trigonal warping. 
Moreover, additional finite energy Dirac points are also generated, but at different 
positions of MBZ in contrast to Park et. alJ 3 ^ Since the SL potential also breaks 
the particle-hole symmetry, the energy spectrum no longer possesses this symmetry, 
and only in an energy window below the original Dirac point, the newly generated 
massless fermions are truly Dirac fermions and the density of states can become 
zero, while for those above, massless fermion states are obscured by the presence of 
other states. Recently, a scanning tunnelling microscope measurement of graphene 
on hBN has observed dips in the differential conductance and thus confirmed the 
existence of finite energy Dirac points.^ 

3. Superlattices in bilayer graphene 

We now turn our attention towards ID electrostatic potential modulations in BLG. 
In general, the features of the band structure will depend on the details of the 
superlattice potential^ 53 * 54 !. For the bilayer system, a general modulation can be 
decomposed into two basic types: i) a chemical potential modulation where both 
layers sit at same potential and, ii) an electric field modulation where there is a 
local interlayer bias. If the SL potential is purely of one type, there is a dramatic 
restructuring of the band structure, particularly at low energy. Notably, the low en- 
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ergy quasiparticles transform from being massive chiral fermions in intrinsic BLG 
to massless chiral Dirac fermions for certain SL parameters. In both cases, much 
of the band structure can be understood by appealing to the inherent symmetries 
and/or to an intuitive effective low energy model. The generation of new zero energy 
modes has similarly been shown to arise in a periodic array of ^-function poten- 
tials^, twisted BLG^, and along domain walls in monolayer graphene with broken 
sublattice symmetry* 56 * 57 !. 

In this section, we focus on reviewing the band structure of the two rudimentary 
types of SL in BLG, a chemical potential and electric field superlattice. Both types of 
SL are of particular interest because each can support the formation of new Dirac 
points for arbitrarily weak SL strengths, in contrast to SL in the monolayer. In 
fact, the Dirac points for the electric field SL survive even for strong modulations. 
A thorough understanding of these two basic SL potentials also provides a firm 
foundation to understand more generic SL profiles and the formalism reviewed here 
can readably be applied to more general SLs. 

We start here by introducing the low energy Hamiltonian that can be used to 
study the properties of generic SL of moderate strength. It should be noted that for 
larger SL potentials, the full tight-binding model is required to correctly describe 
new features in the band structure. Instances where the full Hamiltonian gives 
quantitative differences in the band structure will be duly noted. After establishing 
the formalism, we discuss the band structure generated by a chemical potential and 
electric field SL in Sec. |3.1| and Sec. |3.2[ respectively. 

The low energy Hamiltonian for Bernal-stacked BLG can be obtained by ex- 
panding its minimal tight binding spectrum near one of the Brillouin zone corners 
(K points).^ When the layer potential (i.e., interlayer potential difference) is not 
too large, |A| <C t±, we find H = ip^Hip, ESl where 

fi = _v^f (sp x +ip v ) 2 
tj_ \{sp x - ip y ) 2 

and ip T — ( fl x,frx), with a (b) being the electron operator on the top (bottom) 
layer. Here, p x t y ) — —id x ^ is the momentum operator, ,s = ±1 for the Hamiltonian 
at the ±K valley, vf = V3td/2 10 6 m/s is the Fermi velocity, i«3 eV is the 
nearest neighbor hopping integral, g?«2.46 A is the distance between neighboring 
atoms on the same sublattice (note: d = a\fi where a is the nearest neighbor 
Carbon-Carbon distance), Vi,2 are the potentials on each layer, and t± « 0.15t is 
the interlayer coupling. Unless stated, we set t = d=l. We will ignore inter-valley 
scattering assuming the potentials are varying slowly on the scale of d, so we only 
consider the s = +1 valley (at K). Such an approach has been successfully used to 
study SLs in monolayer graphen d 31 * 32 l 

To diagonalize i/kin, we Fourier transform and then make a unitary transforma- 
tion a p = (a p +/3 p )/v / 2,b p =e 2ie v {a p ~/3 p )/V%, where cos 9 p = p x /p and p = Jp\ +p%. 
This leads to H kin = J2 P ( £ e(p)(3 P f3 p +Sh(p)a p u p ) . Here e eyh (p) =±p 2 /2m* are en- 



)+r ( oV,j. « 
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ergies of electron (hole) states, with an effective mass m* = t±/(2vp). This minimal 
model supports quadratic band touching points at ±K. 

When Vi t 2(x) are periodic, we can also Fourier transform the SL potential to 
obtain H SL = £ p G ^(p)W p ^{p - G), where 

1 fV 1 (G)+V 2 (G)e 2ie V 1 (G)-V 2 (G)e 2i ' 
pG 2 \V 1 {G)-V 2 (G)e 2ie V 1 (G) + V 2 (G)e 2i < 

^(p) = (aj,, (3^), and 9 = 9 p ^g—9 p is the angle between momenta p— G and p. Our 
aim is to understand the band structures of SLs described by -ffkin + ^sl- We will 
study ID SLs with period A along y, so that the reciprocal lattice vectors, {G}, are 
integer multiples of Q = (0, 2tt/X), and the mini Brillouin zone (MBZ) boundaries 
are at p y = ±ir/\. 




3.1. Band structure of ID chemical potential superlattices 

A chemical potential SL corresponds to the case where V\(x 1 y) = V 2 (x, y) = U(x, y). 
For simplicity, we first consider a step-like potential with with (i) U(x,y) = U for 
< y < A/2 and (ii) U(x,y) = —U for A/2 < y < A and use the effective two- 
band Hamiltonian introduced above. Starting from U = and increasing the SL 
strength to moderate values, we observe the following restructuring of the band 



dispersion (see Fig. 11): i) the zero energy quadratic band touching point splits 
into two anisotropic Dirac cones located at (0, ±p y ), ii) further increasing U causes 
the Dirac points to push out towards the MBZ and, iii) upon reach the boundary 
at (0,±7r/A), a band gap opens at a critical U = ?7 e J 28 l 39 l 5 21 Before considering the 
band structure for U beyond U c , let us first discuss the formation of the Dirac cones 
in more detail. 



3.1.1. Dirac Cones: Formation 

The sequence of semimetal to band insulator with increasing SL strength was shown 
to not be dependent on any symmetry in SL profile and to be robust even against 
weak perturbations that vary slowly perpendicular to the principle SL direction.^ 
Given the persistence of the Dirac point, it important to understand why it forms 
and how it is protected. 

The formation of linear band crossing points has been argued to be deeply 
rooted in the chiral nature of the low energy BLG quasiparticles. 1221221 This can be 
seen from the scattering angle dependence of the matrix elements in Eqn. [24] that 
arises from the pseudospin structure of the eigenstates. For states with momenta 
parallel to the modulation direction, 9 = or w, the off-diagonal matrix elements 
vanish; the electron and hole states decouple, so that a particle in an electron (hole) 
state can only forward/back-scattering of the SL potential to another electron (hole) 
state. Since all such electron (hole) states within the first MBZ in an extended zone 
scheme only mix with electron (hole) states of higher (lower) energy, the energy 
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of the conduction (valence) band will be globally shifted down (up). This results 
in two level crossings along the modulation direction, which are protected by the 
chirality of the low energy BLG quasiparticles. If this electron-hole decoupling was 
true for all momenta, we would see the two parabolic bands crossing on a full circle 
in the MBZ, but going to momenta (8p x ,p y ) leads to electron- hole mixing that is 
linear in dp x ; this results in an avoided level crossing and the robust emergence of 
two Dirac cones in the MBZ. 



\ f 

2 14 2.1 2.06 3J58 3.62 3.66 

k k 




Fig. 11. Energy spectrum for a ID superlattice with step-like chemical potential modulation of 
amplitude U. We set A = dOd, with [left panel] U = O.Oli showing two Dirac nodes split along y 
near K, and with [right panel] U = 0.04t showing a full gap. 



3.1.2. Dirac Cones: Properties 

The location and velocity anisotropy of Dirac cones, as well as the critical mod- 
ulation amplitude to gap them out, can be estimated using perturbation the- 
ory in U(G). The second order energy correction of states with p = (0,pj,) is 
A£« 2 »(p) = En^ol (7 («Q)l 2 /[e e ,/ l (p)-e e ,/ l (p + nQ)]- Since e e (p) < e e (p + nQ) 
while £/i(p) > Ship + n Q) m the MBZ, this correction is always negative (positive) 
for electron (hole) states, as expected. 

Thus, the two bands will intersect and cross linearly at momenta (0, iPy), where 
pf/2m* = 2m* £„ #0 \U(nQ)\ 2 / [n 2 Q 2 + 2p* y nQ] . For weak modulations, p* y /Q< 
1, and keeping only n = ±l, we estimate p* « v2to*|£/(Q)|A/7t. For a step profile, 
|E7(Q)| =2U/tt, and H>1 contributions are small. 

For small 8p x away from the level crossing point, we can estimate the electron- 
hole mixing term using perturbation theory and we hnd that the resulting eigen- 
states have energies e p = ±(16m*\U(Q)\ 2 /\Q\ 2 )5p x /p*. The crossing points at 
(0,±p*) are thus really massless Dirac points in the full MBZ. We find velocities 
v y =p*/m* R3 v / 2A|C/(Q)|/7T, and v x = 2v y for the anisotropic linear dispersion. 

Once these Dirac nodes reach the MBZ boundary, Bragg scattering between 
them opens up a full gap. The critical potential strength, |U C (Q)| for this is roughly 
estimated by setting p* = Q/2, which yields |£7 C (Q)| s» tt 2 / (\/2m* \ 2 ) . For a step 
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profile, with A = 60d, we find U c « 0.03i which is close to the numerical result 
0.021 



3.1.3. Strong Potentials and Electron Screening 

For larger SL potentials beyond U c , it becomes necessary to employ the four band 
model to capture a the effects of the high energy bands. Increasing U above U c , the 
band gap continues to grow until a maximum value is reached and then begins to 
decrease before finally closing by forming two new pairs of Dirac cones. For even 
greater U, two of the four Dirac cones merge at (0, 0) and become gapped, but 
remaining two Dirac cones retain the semi-metallicity of the system.'^SMl 

Tan et al.^also performed a self-consistent tight-binding calculation to describe 
the higher energy effects of the entire band structure and to determine the effects 
of interactions on the band dispersion. Besides confirming that the simple single 
particle low energy model correctly describes the qualitative features of the band 
structure, it showed that it is possible to account for screening at the Hartree level 
by a dielectric constant e ~ 11. Hence, the main effect of electron interactions is 
to screen the external SL potential, therefore increasing the critical SL potential 
required to open a band gap. 



3.2. Band structure of ID electric field superlattices 

When BLG is subjected to an electric field SL the potentials on the two layers are 
such that V\[x,y) — —V2(x,y) — U(x,y). In contrast to the chemical potential SL 
discussed above, the band structure is sensitive to the form of the SL profile^. 
To illustrate this, Killi et al."™ considered a more general periodic potential, with 
U(y) — 2U(1 — w/A) for < y < w, and U(y) = —2Uw/X for w < y < A, where we 
have kept the average potential to zero. 

If the parameter w = A/2 the SL potential is symmetric. A numerical calculation 
of the band structure show a pair of anisotropic massless Dirac cones forming at 
zero energy at (±p*,0), as seen in Fig. 12 (left panel)^^. Here, the zero energy 



Dirac cones lie along the direction perpendicular to the modulation as opposed to 
along it for chemical potential SL. Two additional anisotropic Dirac cones are also 
present at high energy, one in the valence band the other in the conduction band. 
Irrespective of the strength of the SL potential, these Dirac points are pinned MBZ 
boundary at (0,7r/A) (or equivalently (0, — ir/\)). 

When w ^ A/2 a band gap opens at all of the Dirac points. This suggests that 
the protection of the Dirac points is governed by a symmetry of the SL profile that 
is broken when w =/= A/2. It is found that the relevant symmetry corresponds to a 
generalized parity operator V that transforms y — > —y followed by exchanging the 
two layers of BLG, and the Dirac points persist as long [P, H) = 0. 

In the letter by the present authors,^ a simple intuitive picture for understand- 
ing all of the features of the Dirac cones was proposed. The idea is to view the 
SL potential as establishing a periodic array of 'kink' and 'anti-kink' steps in the 
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Fig. 12. Energy spectrum for a ID symmetric (see text) electric field superlattice with A = 60d 
and U = 0.03t, showing a pair of zero energy massless Dirac fermions at (±pj,0) [left panel] and 
a nonzero energy Dirac point at (0, ±7r/A) [right panel]. 

potential profile where the parity of interlayer bias reverses. Since it was shown by 
Martin et al.^ that topological zero energy modes are confined along an isolated 
kink (or anti-kink) , it is possible to construct an effective low energy theory evolving 
these modes. Moreover, the band structure of the SL should be entirely dictated by 
these modes fore energies below U where the bulk states are gapped. As a prereq- 
uisite to presenting and analyzing the low energy effective model, it is necessary to 
understand some of the basic properties of the ID kink modes. The next we provide 
a brief overview of these states before returning the construction of the low energy 
model of the electric field SL. 



3.2.1. Kink in the electric field: Soliton modes 

For a general potential profile with V g (y > 0) = — V g (y < 0) and V g (y — > ±oo) = 
±V g , the bulk region far from y = has a gap A « V g , while along this interface, 
localized 'topological' edge modes emerge that are analogous to those in quan- 
tum hall systems.^ These modes can be thought of as forming chiral ID quantum 
wires, since states from opposite valleys are counterpropagating (this follows from 
the Berry curvature about each valley having opposite sign). With respect to the 
two band Hamiltonian, the kink interface marks a region where the mass of the 
quasiparticles, oc a z , changes sign. 

Solving the full tight binding model for a single kink yields the dispersion de- 



picted in Fig. 13 A single kink interface generates right moving subgap modes in 
one valley and two left moving modes in the opposite valley (labeled in red). For an 
anti-kink profile, the dispersion is identical except the modes velocities are reversed 
in each valley. 



In terms of the low energy, the eigenfunctions are then of the form 

(f(y)\ = ( m \ ( f(v) 



\g(y)J 0/n \f(-y)J \-f(-y) 



(25) 



with corresponding eigenvalues of +1 and —1 of the operator V, respectively. So- 
lutions with eigenvalues +1 with even symmetry belong to the lower 0-band while 
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solution with odd symmetry belong to the upper 7r-band. 

Increasing the interlayer bias strength results in two important effects that can 
be seen qualitatively in Fig. [l3j i) the Fermi velocity of the two bands is enhanced, 
and ii) the wavefunctions become more confined to the interface. With respect to 
the overall width of the wavefunction, a simple scaling analysis suggests that the 

wavefunction width should go as I ~ (m*V g )~ 1 ~ I -A — 




- -\ V n 



4.3 -4jt/3-4.I ', ' 4.1 4n/3 4.3 -4.3 -4jt/3-4.l ', ' 4.1 4lt/3 43 -100 









_ v e =0.02t 








_ V =0.08t 


l())l 2 






1 




W 




\ 2 



100 



Fig. 13. Dispersion about the K-points with (a) V g = 0.02t and (b) V g = 0.08t. Edge-mode bands 
are indicated by the labelled arrows and bulk-states by the hatched region. The modulus square 
of the zero-energy wavefunction of (c) the 0-band at V g = 0.024 and V g = 0.08t (the 0,2 and b\ 
branches are exchanged for the 7r-band). 



The large width of the wavefunction transverse to the wire direction strongly 
suppresses the bare backscattering terms due to electron-electron interactions. A 
similar effect also seen in wide carbon nanotubes,^ and it leads to a dominance 
of forward scattering processes, where a simple bosonization analysis predicts a 
spin-charge separated gapless Tomonaga-Luttinger liquid. We thus expect a rela- 
tively large energy window where interactions drive the ID modes even in these 
'kink' modes in bilayer graphene into such a Luttinger liquid. Remarkably, such 
a bosonization analysis arrives at a novel two-band Luttinger liquid with tunable 
mode velocities and tunable Luttinger parameters.^ 

Naively, it would appear that the localized kink states may not be robust be- 
cause they are not topologically protected. However, Qiao et al,^ performed an 
extensive study into various potentially detrimental mechanisms. Through numeri- 
cal conductance calculations and examining the LDOS, they showed the low energy 
states are remarkably robust to both short and long range disorder, and even to 
abrupt changes in the interface direction. Although quantized conductance is not 
unlikely, the mean free path could be as large as a hundred microns for relatively 
clean samples. The mechanism which quantitatively leads to a strong suppression 
of backscattering is again the large wavefunction spread. 

It has also been observed that the kink-modes are also quite robust when sub- 
jected to a magnetic field.^^ This is due to the strong magnetic field induced 
confinement of the wavefunctions J^Sl Interestingly, by coupling a pair of coupled 
kink and antikink modes and applying a magnetic field, the current in the kink 
flows in one direction, opposite to the direction of flow in the antikink. Moreover, 
all the dispersing modes have the same valley index, making this a potential valley 
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filter .1441SQ1 Recently, it has been proposed that at v = electron interactions form a 
charge density pattern in the vicinity of a kink state, which provides a key signature 
of quantum Hall ferromagnetism.SS 

Before turning our discussion back over to superlattices, we close this section by 
emphasizing that similar localized kink states are also expected to form naturally 
in the presence of charge impurities and in strongly correlated phases where the Z2 
layer symmetry is spontaneously broken. In the case of the former, charge impurities 
close to the surface of the sample can generate a local electric field strong enough 
to induce an interlayer bias. In uniformly biased BLG or where there are multiple 
charge impurities in close vicinity (on opposite layers or with opposite charge), 
this can cause the interlayer bias to reverse, generating kink states - a point to be 
elaborated on in the conclusion. In the latter case, any state with spontaneously 
broken layer symmetry will naturally form domain walls separating regions with 
opposite interlayer bias .121121 Again, kinks states are expected to from percolation 
networks that permeate throughout the bulk. 



3.2.2. Electric Field Superlattice: Effective Model 

Equipped with an understanding of the properties of the soliton kink/anti-kink 
modes, it is now possible to describe how to construct a low energy effective model.^1 
To begin, first consider the dilute limit where the period length is much longer than 
the characteristic spread of the soliton modes, (i.e. A >> /). Each kink supports two 
(ignoring spin) unidirectional dispersing soliton modes while each anti-kink supports 
two oppositely moving modes per valley, as shown in Fig. 14 The counterpropagation 
of the kink and anti-kink modes results in four band crossing points about each K- 
point, two at zero energy between bands with the same symmetry (0-0 and tt-tt) 
and two at finite energy between modes with opposite symmetry (w-0 and 0-7f). As 
described below, when the wavefunctions of neighbouring soliton modes couple (i.e. 
A I), Dirac cones precipitate precisely at the band crossing points. 

At energies and momenta in the the vicinity of any one of the band crossing 
points the system looks as if it were an array of ID chiral 'wires' lying along the 
kinks and anti-kinks of the SL. Each wire supports modes that flow in opposite 
direction to its two neighbors. Now, as the wavefunctions of these modes begin to 
overlap the electrons can hop between neighbouring wires. 

With this in mind, let us consider the ir-0 modes at zero energy and at a mo- 
mentum p* (away from K). The hopping between neighboring wires along y is then 
between states which have opposite velocities (since it is between a kink and an 
anti-kink edge state) and it is between a p-wave like state (P-odd) and an s-wave 
like state (P-even) (see Sec. 3.2.1). 

Careful attention must be made to get the correct form of transfer integral that 
describes the hopping between the wires. The sign of the hopping can be deduced 
by taking the wave functions of the 0-band and 7r-band as having s-wave orbital and 
p-wave orbital character, respectively, and noting the sign of the overlap between 
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the wires. An illustrative picture for two of the band crossing points is shown on the 



right in Fig. 13 Hence, the hopping between wires at a zero energy band crossing 
point is staggered and uniform for the finite energy band crossing points. Further 
details are provided in our previous work.^ 

Using the index n to label the wires, the interchain hopping parameter will then 
alternate as (— l) n g for equally spaced wires and as g + 6, —g + 8 (with S < g) if 
pairs of wires are closer to each other. Linearizing the dispersion at the crossing 
point, and letting vq denote the velocity of the linearized modes, 

H(p x ) = v -Pl)4*n c P* n ) 



{g(-l) n + S) (4.n9p.n+l + h.c.) 



(26) 



where p* is the location of the n — crossing point in the single kink or antikink 
problem, and c Px „ annihilates an electron on wire n with momentum p x . Let t;(p x ) = 
vq(p — P* x )- Fourier transforming, we find H{ Px ) = Yl Py &(p y )<T-h(p x )*(p v ), where 
HPx) = {£{Px), -2gsin(p y ), -2Scos(p y )), with *(p y ) = (c Py c Py+7r ) T , and Y! Py runs 

over the MBZ. The dispersion is thus E = ±^J( 2 (p x ) + 45 2 cos 2 (p y ) + Ag 2 sin 2 (p y ) . 
Consequently, when w = A/2, and the Hamiltonian commutes with V, we have 
5 = and a Dirac cone is generated at (p*,0), consistent with numerical results. 
When w ^ A/2, the Hamiltonian breaks V — we then have 6^0, which leads to 
a gap 4(5. Similar arguments hold for the other zero energy band crossing points. 
The velocity of the Dirac fermions is highly anisotropic and depends on <?, except 
along the SL direction where it inherits its value from the freestanding zero mode 
velocity — this can be controlled by tuning the SL period and amplitude. 

There are a number of particularly salient properties of the electric field SL 
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Fig. 14. (color online) Left: Spectrum of isolated kink (thin, red) and anti-kink (thick, blue). 
Higher (lower) energy modes are labelled it (0) at a kink and as # (0) at an anti-kink. Right: 
Schematic of hopping between the it — n and — tt states. 
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when viewed from the 'coupled wire' perspective. Although for different reasons 
than the chemical potential SL, its band structure is also expected to be quite 
resilient to disorder. As discussed previously, the underlying soliton modes of the SL 
are exceptionally robust to disorder and fortifies the the band structure. Moreover, 
since the individual ID soliton modes are present over a wide range of interlayer 
biases, the Dirac spectrum of SL is persists for both weak and strong SL potentials. 
In addition, the properties of the low energy Dirac fermions are very versatile. 
Specifically, the velocity parallel and perpendicular to the modulation direction 
can be tuned independently by first adjusting the SL strength and then the period 
length. Furthermore, mass can be imparted to the fermions by breaking the V- 
symmetry of the SL. Interestingly, just as in polyacetylene, a domain wall between 
a gapped region with w > A/2 and a gapped region with w < A/2 leads to new 
subgap soliton modes. Since each kink/anti-kink is itself like a domain wall, these 
should be viewed as solitons in a soliton lattice. 



3.3. Magnetic Field Effects on ID superlattices 

In this section, we review the effects of a magnetic field on the single particle prop- 
erties of BLG subject to ID SLs, but before doing so, it is useful to briefly review 
the LLs of intrinsic BLG. The low energy model that describes BLG in the presence 
of a perpendicular magnetic field is obtained by replacing the momentum operator 



in Eqn. 23 with its canonical counterpart to take into account of the magnetic field 
effect. 

The following results were obtained from effective two band model with the 
same gauge choice as before, A — Byx. The eigenvalues and the corresponding 
eigenvectors of the above Hamiltonian for the s = +1 valley in the absence of a SL 
are 



£„ = sgn(n)y/ \n\(\n\ - l)uJ 2 Jt ± , 

V2L V- s g n ( n )V'|n|-2,fe(y)/ 

with \n\ > 2. In addition, there are two zero energy solutions that are feature the 
hallmark feature of the quadratic band touch point OsSED 

£l = o, ^A^-^^f 

e = 0, ^A^^^^f)- (28) 

For s = —1, the corresponding eigenvectors are given by 4> n .k.-{ x iy) — 
<J x't > n,k,A x i V)- The full low energy LL wavefunctions thus take the form 
4>n,k,±e ±iKxX ■ and these serve as a good basis to study the magnetic field effect 
of bilayer graphene SLs. 
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3.3.1. Chemical Potential Superlattice: Landau Levels and DC Conductivity 

Just as for the Dirac cones derived from the SL in the single layer, evidence of Dirac 
fermion dispersion was shown to exist in the LL structure of ID SL in BLG.^The 
LL energy spectrum generated by a weak magnetic field is shown in the left panel 
of Fig. [15] for Vo = O.Oli. Degenerate pairs of energy levels that derive from the 
anistropic Dirac cones are present at zero energy (n = 0) and at finite energy 
(n = ±1). The higher LL, however, do come in degenerate pairs because it is only 
at lower energies the spectrum consists of two copies of Dirac cones. The right panel 
shows that for Vo = 0.04i, the zero energy LL levels are absent, consistent with the 
opening of a band gap. 




Fig. 15. (Color online) Left and center panel: Energy spectrum of BLG subject to a chemical 
potential SL and a weak perpendicular magnetic field (£g = 2A) for Vo = O.Oli and Vo = 0.044, 
respectively. Right panel: Evolution of low lying energy levels as a function of SL potential strength 
U, with lg = 2A, yo = 0. In all cases, A = 100a, where a = 1.42A. 



The evolution of the LLs as the SL potential is increased, shown in Fig. 15 



displays a definitive crossover from a non-relativistic to a relativistic regime, in 
addition to the opening of a bandgap. As SL potential increases, the physics grad- 
ually becomes dominated by anisotropic Dirac cones, which can be seen from the 
appearance of doubly degenerate levels at nonzero energies. The existence of a 
marked crossover can be qualitatively understood by considering the competition 
between the characteristic energy scales in these two regimes. In the absence of 
the SL, the low energy excitations are massive electrons with an effective mass 
m* = t±/2vp and have a cyclotron frequency, ui' c = eB/m*c — l/m*i 2 B . On the 
other hand, the anisotropic Dirac points generated by SL have anisotropic Fermi 
velocities v y = v / 2A|[/(Q)|/7r and v x = 2v y , where Q = y2n/X^ and have the 
characteristic energy scale of uj c — ^/2v x v v /Ib- Wu et al.^ estimated the crossover 
should occur around U ~ 0.002i, which is quite close to the value observed in Fig. [15] 
(right figure). 

Further increasing the SL potential, the doubly degenerate zero energy levels 
become gapped and all levels are pushed away from Dirac point. Surprisingly, at 
rather strong SL potential, U ~ 0.22t, zero energy LLs appear again, and all of the 
higher energy levels become doubly degenerate. This phenomenon can be under- 
stood from the result of Tan et al.^ As it has been shown, for a chemical potential 
SL, when SL potential is strong enough, anisotropic Dirac cones will show up again 
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-So 







Fig. 16. (Color online) Diagonal dc conductivities for bilayer graphene for a chemical potential 
SL with different strengths Vb, and magnetic fields B. The conductivity is shown for weak field 
(£g = 2A, top panels) and intermediate field (£g = 0.2A, bottom panels). Left panels (a,b) 
correspond to Vb = O.Olt, and right panels (c,d) correspond to Vb = 0.044. The conductivities 
show anisotropy, where a xx is always larger than o- yy , in contrast to anisotropy reversal in MLG 
SLs. 



in the energy spectrum, which naturally leads to the zero energy LL in the presence 
of a magnetic field. As shown elsewhere,^ there can be up to four Dirac points in 
the spectrum. For even stronger chemical potential SLs, the degeneracy of the zero 
energy LL reduces to two, consistent with two of the Dirac cones becoming gapped 
as discussed in Sec. 3.2 Fig. 16 shows the dc diagonal conductivity of chemical 
potential SLs, where we again see anisotropy similar to the single layer case. How- 
ever, in contrast to the monolayer, the direction with largest conductivity does not 
reverse when the magnetic field strength is tuned, and so the anisotropy cannot 
be tuned. For weak fields, the transport anisotropy is directly determined by the 
anisotropy of the Dirac cones. In Ref.^, the extent anisotropy was determined by 
treating the SL potential as a perturbation. It was calculated to be v x ~ 2v y for the 
emergent Dirac cones, consistent with the observation that the conductivity in the 
x direction, a xx , is larger than a yy in a weak magnetic field. As for intermediate 
magnetic field, due to the dispersion of the energy bands, the average velocity in 
the x direction is not zero, (v x ) ^ 0. On the other hand, (v y ) is always equal to 
zero. This means that a xx will acquire intra-LL contributions, while a yy is mainly 
determined by inter-LL contributions and is thus small compared to a xx . 
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Fig. 17. Left: Energy spectrum of an electric field superlattice, in a weak magnetic field. Centre: 
Evolution of low lying energy levels in an electric field SL as a function of SL potential strength 
U with y = 0. Right: DC conductivity. In all plots, l B = 2A, , A = 100a where a = 1.42A. 



3.3.2. Electrical Field Superlattice: Landau Levels and DC Conductivity 

In a ID symmetric electric field SL, there are always two zero energy Dirac points 
present in the spectrum, which results from the coupling of ID zero modes at 
kink/antikink of the SL potential profile. Wu et al.^ argued that this implies that 
when a weak magnetic held is applied, doubly degenerate zero energy levels should 
appear at the Dirac point. Indeed, as can be verified from the LLs shown in the left 
panel and the evolution of the LLs in the centre panel of Fig. [TTJ these two levels 
are always present at zero energy and are independent of the SL potential strength. 

The authors also showed that for a SL potential chosen to be Vq = 0.03i, nearly 
degenerate levels even appear at nonzero energies (Fig. 17 1. These correspond to the 
LLs derived from anisotropic Dirac points, up to n — ±4. From the left panel of Fig. 



17 it is more clear that at strong SL potential, physics is strongly dominated by the 



Dirac points, where higher energy levels become doubly degenerate and resemble 
the higher LLs of the Dirac cones. When SL potential is weak, equally spaced LLs 
are recovered, as in the chemical potential SL, which also indicates a nonrelativistic 
to rclativistic crossover at certain SL strength. Remarkably, and different from the 
chemical potential SL case, the relativistic behavior survives to higher energies as 
SL potential increases, which means the linear approximation description of Dirac 
cones works in a larger energy range. This is consistent with earlier result J2SM] 
Therefore, as the SL potential increases, the energy range where the Dirac cone 
approximation is valid also increases, which leads to the robust relativistic physics 
at large SL potential. In addition, evidence for the large anisotropy of Dirac cones 
discussed previously can be readily seen in the DC conductivity shown in the right 
panel of Fig. [17| 



4. Concluding remarks 

We have explored, here, the rich physics associated with slowly modulated potentials 
in monolayer and bilayer graphene. The resulting dispersions and magnetotrans- 
port properties of such superlattices can be explored experimentally by engineering 
gates to pattern suitable superlattice potentials. Indeed, such periodic modulated 
potentials for Dirac fermions have also been explored in recent scanning tunneling 
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spectroscopy studies the surface states of Bi2Te3, a topological insulator, where 
the modulation arises from the periodic buckling of the crystal structured Going 
beyond such slowly varying potentials, there has been a lot of recent interest in 
graphene on hexagonal Boron Nitride substrates, where lattice mismatch leads to 
Moire patterns - such potentials imparted by the hBN substrate have slow as well 
as fast sublattice scale components, leading to new physics which is still being ex- 
plored. EUHMMl Experimentally, tunneling studies indicate the emergence of new 
finite energy Dirac points in this case,^ and further studies in this area would be 
valuable. Turning to a different aspect of such potential modulations, we note that 
the issue of transport in bilayer graphene under a uniform bias is likely to be im- 
pacted by the presence of charged impurities in the substrate. E3HSI1 ^he electric field 
of such an impurity could locally reverse the applied bias, leading to a ring around 
the impurity site which can trap midgap states.^ While the ring size is likely to 
be small, on the nanometer scale, the trapped states would have a large wavefunc- 
tion spread due to the very small effective mass associated with the quadratic band 
touching point in bilayer graphene. In this case, we expect the low temperature 
transport could be via hopping conduction between such "ring" localized states, 
while the observed small activation in biased samples^^l could result from such 
hopping conduction getting gapped out due to the charging energy of such "ring" 
states. Rough estimates show that the wavefunction spread could be on the order of 
50nm, while the charging energy could be ~ lOmeV, comparable to gaps observed in 
transport measurements .IS4IS2I The study of a network model of such "ring" states 
is likely to yield new insights into transport mechanisms in bilayer graphene. In 
summary, the study of periodic and random potential modulations in graphene and 
bilayer graphene is a rich and growing field with many open questions and we invite 
the reader to join us on this exciting journey. 
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